rm(list=ls())

load("TSCS_data.RData")

# Table 3 -----------------------------------------------------------------

lm_support <- plm(citizen_support ~ 
                    lag(citizen_support,1:2)  + 
                    antipluralism_gov_seat_share_lag + 
                    lag(polarization,1) +
                    lag(clientelism_gov_seat_share,1) +
                    lag(region_citizen_support,1) + 
                    lag(region_antipluralism_gov_seat_share,1) + 
                    lag(region_liberal_democracy,1) + 
                    lag(democratic_stock,1) + 
                    lag(GDP_capita,1) + 
                    lag(life_expectancy,1),
                  data = panel,
                  model = "within",
                  na.action = na.omit,
                  effect = "twoways")

cmodel_support <- pgmm(citizen_support ~ 
                         lag(citizen_support,1:2)  + 
                         antipluralism_gov_seat_share_lag + 
                         lag(polarization,1) +
                         lag(clientelism_gov_seat_share,1) +
                         lag(region_citizen_support,1) + 
                         lag(region_antipluralism_gov_seat_share,1) + 
                         lag(region_liberal_democracy,1) + 
                         lag(democratic_stock,1) + 
                         lag(GDP_capita,1) + 
                         lag(life_expectancy,1) |
                         plm::lag(citizen_support, 3:5),
    data = panel,
    effect = "individual",
    model = "onestep",
    transformation = "ld",
    indexes = c("country", "year")
  )

# Table model
modslin <- list(lm_support)
vcm_lm <- lapply(modslin, function(x) plm::vcovBK(x, cluster=c("group")))
rse_lm <- lapply(vcm_lm, function(x) sqrt(diag(x)))
pval_lm <- lapply(1:length(modslin), function(i) {
  pt(abs(coef(modslin[[i]]) / rse_lm[[i]]), df = 1608, lower.tail = FALSE) * 2
})
model1_wool <- sapply(modslin, function(x) round(pbgtest(x, order=1)$p.value, 3))[c(1,0)] # Wooldridge test for serial correlation

modsmom <- list(cmodel_support)
vcm_gmm <- lapply(modsmom, function(x) plm::vcovHC(x, cluster=c("group")))
rse_gmm <- lapply(vcm_gmm, function(x) sqrt(diag(x)))
pval_gmm <- lapply(1:length(modsmom), function(i) {
  pt(abs(coef(modsmom[[i]]) / rse_gmm[[i]]), df = 1608, lower.tail = FALSE) * 2
})
model1_AB <- sapply(modsmom, function(x) round(mtest(x, order=2, vcov=plm::vcovHC(x, cluster="group"))[["p.value"]][[1]] ,3))[c(1,0)] # AR(2) p-val
model1_H <- sapply(modsmom, function(x) round(sargan(x, weights="twosteps")[["p.value"]][[1]], 3))[c(1,0)] # Hansen test p-value

stargazer(modslin,modsmom,
          covariate.labels=c("Citizen Support t-1","Citizen Support t-2",
                             "Anti-pl. government t-1","Polarization t-1","Clientelism t-1",
                             "Region Citizen Support t-1","Region Anti-pl. Government t-1",
                             "Region Democracy t-1","Democratic Stock t-1","log(GDP t-1)","Life Expectancy t-1"),
          type="latex",
          omit.stat	= c("f","n"),
          se = list(rse_lm[[1]],rse_gmm[[1]]),
          p = list(pval_lm[[1]],pval_gmm[[1]]),
          dep.var.caption	= "",
          dep.var.labels	= "Citizen Support",
          add.lines=list(c("N Observations", "1895","3691"),
                         c("N Countries", "99", "100"),
                         c("Wooldridge test (p-value)", model1_wool,""),
                         c("Arellano-Bond test (p-value) ", "",model1_AB),
                         c("Hansen test (p-value)", "",model1_H)
          ),
          no.space = T,
          star.cutoffs = c(0.1,0.05,0.01),
          out = "tab_3.tex"
)